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1 Introduction 

We propose and study a method for a multi-agent system of au- 
tonomous vehicles to perform the exploration of a non-stationary 
environment. The goal of this exploration is to reach a desig- 
nated target geometry, typically a boundary, and to describe that 
boundary. The current paper relates to the general objective of 
studying motions of swarms, which we understand as a collection 
of autonomous entities which rely on local sensing and simple 
behavior, interacting in a way that a more complex behavior of 
the whole group emerges from local interactions jSj- This type 
of behavior is well-known in natural phenomena such as schools 
of fish [Sj and colonies of insects One can find in the litera- 
ture various other methods of simulating motions of distributed 
systems inspired by biology |1] , [Sj , and by economics-based con- 
cepts |j7j . The applications of such motions are quite broad since 
substituting agents for humans is desirable in many areas such as 
de-mining operations jBj, mapping and exploration |H], navigation 
control [9j , formation flying fn]\ and military planning . 

The algorithm presented here models a group of agents to move 
on a given surface or within a prescribed volume, find its bound- 
ary and track it, while communicating with each other and spread- 
ing simultaneously over the surface. The environment in which 



the motion takes place varies in time. Such a method may pro- 
vide insight to a better understanding of motions of swarms in 
biology [i4j. The motion of swarms in military applications has 
been considered in recent years [3j, since it offers an unmanned, 
dependable, flexible, self-organized type of network to achieve 
various military objectives, while replacing the burden of com- 
prehensive communication, since the swarms interact primarily 
with their neighbors. 

The particular goal of identifying and tracking the boundary of 
a surface which we address here, is analogous to the problem 
encountered in image segmentation of tracing the boundary of 
an object. Tracking contours as a function of time has been ad- 
dressed in the pioneer work of Kass pjjj , in relation to localizing 
features in an image. Here we add more complexity to the situ- 
ation in [TT] since the motion of the agents starts far away from 
the boundary. Moreover, once the boundary is located the con- 
tour of agents moves along that boundary rather than keeping a 
fixed location. 

The idea of adapting image processing techniques to simulate 
motions of swarms was proposed in [i2jj , where the image segmen- 
tation technique of locating a contour with a model known as a 
snake [m, or energy-minimizing curve, is adopted. We ex- 
tent that analogy by allowing the formation of agents to move on 
a time-dependent surface or on a surface whose boundary travels 
in time. In the current work we explore the problem of simulat- 
ing the motion of agents which effectively locate and follow the 
boundary of a surface of interest, assuming that the boundary is 
known to each agent and it is moving in time possibly along with 
the entire surface. The motion of agents takes place in either 2D 
or 3D. 

The idea of the algorithm we present is to model the group of 
agents as a contour that deforms towards the boundary of the ob- 
ject. Interactions between nearby agents are assumed, the model- 
ing of these interactions being analogous to having forces acting 



2 



between nearby particles as in an elastic band. An energy func- 
tional results from this physical model in a continuum limit. The 
desired motion of the agents is obtained by minimizing this func- 
tional to which an appropriate goal dynamics term is added, so 
that its minimum is obtained as the agents land on the proposed 
target. 

The algorithm assumes that each agent has the position of the 
boundary, has information on nearby agents, and evaluates its po- 
sition on the surface as well as the local gradient of that surface. 
The movement of each agent can be modeled by a partial dif- 
ferential equation for the velocity of each vehicle. By discretizing 
this equation we derive a numerical scheme for the motion of each 
agent. Each equation contains a term accounting for movement 
along the steepest descent on the surface and a term correspond- 
ing to movement parallel to the tangent to the boundary that is 
being tracked. Coordinated behavior is created by imposing lo- 
cal interaction rules among agents. These terms also dictate the 
manner in which the formation will converge towards the desired 
goal. 

The paper is organized as follows: In the next section we present 
an algorithm for an individual agent to accomplish the desired 
task. Then we continue in section III with adding a sparsing term 
which connects the motions of nearby swarmers and prevents col- 
lisions of agents in real applications. Obtaining collective motion 
local communication rules is essential. Typically these rules de- 
termine a general pattern of significance in applications j3] . Con- 
sequently in section III we present a PDE model that includes in- 
teraction between nearby swarmers. These interactions are based 
on modeling the swarm as the discretization of an elastic string 
with elastic and bending forces acting between nearby particles. 
In section IV we give a numerical method based on the model 
in the previous section. We continue in section V with the cor- 
responding algorithm in 3D, which means the motion now takes 
place on a surface defined in 3D rather then 2D. We end with a 
conclusion section. 
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2 The Collective Motion Algorithm 



The simplest version of the algorithm consists in designing mo- 
tion terms that would lead a single vehicle to move towards a 
given boundary and then move along that boundary. The compo- 
nent of the motion that determines the approach to a prescribed 
boundary consists in moving along the steepest descent direction 
along a surface. To make the vehicle describe the boundary we 
include a term that is parallel to the tangent to the boundary 
of the surface. The velocity of a vehicle is thus calculated from 
these two components giving a system of two ordinary differential 
equations. 

Let C{x, y) denote the concentration function of the environment 
at the agent's position {x^y). By "concentration", we mean a 
function which describes any spatially varying variable, such as 
density, temperature, etc. Let P{x,y) = f{C{x,y)) = (C — Cq)^ 
be a function that achieves a minimum at the boundary of the 
environmental concentration. Let v = v{x,y) denote the position 
vector. The motion according to the rule ^ = —VP results in 
a gradient descent toward a local minimum of the function P, 
therefore towards the curve C = Cq. Also the rule ^ = V^C, 
(where V"^C denotes the normal vector to the boundary of the 
surface) results in travel along the boundary of the concentration 
surface C. Combining the two terms mentioned above we arrive 
at the following model: 



dx 
'dt 



= -aP 



IVCI 



dt 



(1) 



The parameter cu in Eq. [T] determines the speed of the vehicle in 
the direction of the tangent to the boundary. Following [2J as a 
test case, we consider as the concentration function: 



4 



Fig. 1. The trajectory of a single agent (represented by '*'s) launched on a concen- 
tration surface, at initial condition xO = 0.3 and yO = 0.5, reaching the boundary 
and accurately following the boundary. Solid lines denote the surface C{x,y). 
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Equation Elis defined on a disk of radius 2 and {xi^yi) = (1,0), 
(0,-|),(-|,i),(^, 1) for I = 1,2,3,4 and a = 1.0. Co is the 
concentration on the boundary, so it is obtained from Eq. ( [2] ) 
taking x and y on the boundary of the circle of radius 2.0 centered 
at the origin. In this way we are targeting a boundary that varies 
in space. 

We solve the Eq. ( [T] ) subject to Eq. E] by using Runge-Kutta 
method in time to obtain the trajectory shown in Fig. [H where 
the initial conditions were xq = 0.3 and yQ = 0.5. We also show a 
top-down view of the same trajectory in Fig. [2] where the level 
set of the surface are shown as well as the boundary of the surface 
(the circle of radius two) on which the agent moves. 

To achieve movement along the boundary we calculate Co using 
Eq. ( El ) for x and y corresponding to coordinates on the circle 
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Fig. 2. Top-down view of the surface and trajectory in Fig.l. Level sets of the 
projected surface are shown. The boundary of the surface is projected onto the 
circle of radius 2. 

of radius 2 where the surface is defined (any other curve in x — y 
space could have been considered). 

The algorithm can be further extended to allow a vehicle to tar- 
get a boundary which moves in time. In this case, Eq. (ID) is used 
again, the only difference being that Cq is now a function that 
varies in time as well in space. In Fig. El we show three nearby 
trajectories following a boundary which moves in time at a con- 
stant speed vq = 2. Four different instances in time are shown at 
t = 0.6, 1, 1.5 and t = 2.0. 

The existence and uniqueness theorem for ordinary differential 
equations |T3| guarantees that the trajectories of this motion will 
not intersect. In realistic situations when agents are used, collision 
among agents can become a concern. For such a situation we add 
a sparsing term in Eq. [T] meant to repel nearby vehicles fHj. This 
term varies the velocity of a vehicle according to the following 
equations: 




dvi 
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t=0..2 



t=1.0 




Fig. 3. Three nearby trajectories evolving on a concentration surface whose bound- 
ary expands in time at a constant speed. The four figures show the trajectories at 
t = 0.2, 1.0, 2.0, 3.0. Sohd hues denote the evolving concentration surface. 



where the repulsion kernel is given by U{v,w) = Crexp{— 
Here the parameter Ir is a length scale. In what follows Cr = 10, 
and Ir = 0.5. This term is set non-zero only when the distance 
between vehicles is less than a given cut-off distance R; i.e., a 
maximum distance at which the sparsing is active. In Fig. HId) we 
show 4 trajectories whose motion is given by the Eq. [T]to which 
the sparsing term in Eq. El has been added. For comparison the 
same trajectories without the sparsing term are shown in Fig. 



7 




Fig. 4. Four trajectories, on the concentration surface starting with the same initial 
conditions, (a) There is no communication between agents. The trajectories evolve 
in time independently, (b) Similar to (a) but with a sparsing term added into the 
equations of motion, so that agents that get too close repel each other. 

H^). When sparsing is added, we see the 4 trajectories, shown 
in different colors, intermingle several times on the surface until 
smooth trajectories are established. This differs from the orderly 
trajectories that can be seen in Fig. El 

The above algorithm is used when only a few independent vehicles 
are available. When a sufficient number of vehicles are available, 
direct interaction among nearby vehicles can be introduced. This 
accounts for communication among vehicles in real situations and 
for creating desirable patterns of the motion of a swarm, since a 
swarm motion is essentially defined as creating various overall 
patterns and behaviors out of local interaction rules. In what 
follows the approach will be to view the group of agents as a con- 
tinuous curve moving and expanding in time. A nice continuum 
mechanical analogy occurs when the motion is being thought of 
as an elastic band that expands and bends aiming to fit around 
a given boundary |TT]. This approach will be shown in the next 
section. 
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3 Modeling the swarm behavior as an energy minimizing curve 

The algorithm we develop is meant to be used by vehicles with 
sensing capabilities moving through a medium of variable con- 
centration and aimed at detecting a certain object such as the 
boundary of a surface. Interaction between vehicles can be con- 
ceived in various ways leading to a variety of behaviors of the 
overall swarm. In what follows we adopt the point of view of con- 
sidering the group of agents modeled as an elastic object with 
elastic and bending forces acting between nearby agents 
Achieving the desired goal amounts to solving a minimization 
problem. We summarize the related theory of this approach in 
this section and give the corresponding numerical discretization 
in the next section along with numerical examples. 

The collective motion is modeled as the motion of a deformable 
contour which expands and wraps around the boundary, then 
moves around the boundary. Consequently it will be represented 
by a mapping: 

n = ^ R^,s^ v{s) = (x(s), y{s)) 

where is the periodic unit interval. Periodic boundary condi- 
tions are assumed since any agents in the swarm are considered 
to be identical. The algorithm consists in minimizing a certain 
functional E over a space of admissible deformations F. This 
functional E : F ^ R represents the energy of the contour and 
has the following form: 

E{v) = J^[w,\v,{s)\^ + W2Ms)\'^ + P{v{s))]ds, 

where the subscripts denote differentiation with respect to the 
Lagrangian parameter s which is the arc length, and P is a func- 
tion of the environmental data. P also is designed to incorporate 
the goal dynamics we are interested in; i.e. the boundary of a 
surface. The other two terms of the functional E account for the 
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mechanical properties of the contour of agents; i.e., elasticity and 
bending. The first term makes the contour have stretch and the 
second term makes the contour bend [T^. The choice of the pa- 
rameters wi and W2 determine the elasticity and the rigidity of 
the contour. 

If -u is a local extremum of E it satisfies the associated Euler- 
Lagrange equation: 

d,{wiVs) - dss{w2Vss) + VP = (4) 

with periodic boundary conditions. As before let P = {C — Cq)^. 
The relationship between potential functions and desired goals is 
discussed comprehensively in [12] and [TT] . 

Since the collective motion we are designing takes place in time, 
we approximate the solution of Eq. ( HI) as the steady state of the 
following partial differential equation for the contour v. 

dtv = ds{wiVs) - dss{w2Vss) + VP, (5) 

in which the right hand side is the Euler-Lagrange Eq. HI This 
equation describes how each point on the active contour v should 
move in order to minimize the functional E. In the next section we 
show how this equation can be discretized and used in numerical 
examples. 

4 Collective Motion Algorithm along a Virtual Contour 

In this section we present an algorithm for collective motion based 
on the minimization approach presented above. This algorithm 
is obtained discretizing Eq. ( [Sj) and considering each spatial dis- 
cretization point as representing an agent. Thus a group of agents 
will evolve according to a time-dependent scheme and will move 
along a virtual contour which is an approximation of the curve v 
solving Eq. ( EI) . 
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We rewrite Eq. ( 0) as follows: 



d^x oix gg jSx gggg -\- Oj;P 

(6) 

dty = aygg - Pygggg + dyP 



where as previously P{x,y) = \{C{x,y) — Co) p. Using finite- 
differencing to discretize the derivatives with respect to the La- 
grangian parameter s, we obtain: 



a (5 

(7) 

dtVi = -^{Vi+i - + yi-i) - j^{yt-2 - + 6y^ - %+i + y,+2) 

+dyP{xi,y^) 



This discretization gives the equation of motion for a group of 
agents occupying the positions (xj, yi). We define thus the motion 
of the mobile agents by a set of coupled ODEs for the positions of 
the agents. In what follows we take a = 0.01 and P = 0.0001 
as in [2]. 

The equations of motion given by ( [Tj) will evolve agents towards 
the boundary of the concentration surface. Once on the boundary, 
movement along the boundary is generated by adding a term 
parallel to the normal to the surface as in ( [T]). This results in 
the following scheme: 
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a (5 

dtXi = - 2Xi + Xi^i) - J^{Xt~2 - 4x^-1 + 6Xi - 4:Xi+i + X^+2] 



(8) 



a f3 
+dyP{x^, yi) + j"^/^''^'^ + E fe, Vj))- 



The last term represents the sparsing term that has been added 
as in Eq. El In the examples below we let a; = 2.0 and the time 
derivative in Eq. El has been discretized using the forward Eu- 
ler method. With this time discretization it follows that one and 
respectively two-nearest neighbor communication will occur, as 
seen from the discretization of the diffusion term and the dis- 
cretization of the 4th order derivative. Stability issues as well an 
implicit version concerning this algorithm have been discussed in 
[T| where it is shown that stability of the forward Euler scheme 
requires the condition At < g^^^a/i^ ? where h stands for the dis- 
tance between nearby agents. 

Using Eq. ( jHj) one can effectively target and track a level curve 
of a concentration surface. We have allowed in Eq. ( |HI) the con- 
centration on the boundary to be a function of x and y. In Fig. 
(SI we show a swarm of 25 vehicles at four different instances in 
time, each one 20 time steps apart with a time step At = 0.01. 
We show the swarm at the beginning of the motion which starts 
in the center of the surface, then we show three later pictures of 
the same swarm as it circles the surface boundary. Fig. (HI is a 
top-down view of Fig. in which the projection of the surface 
boundary, (i.e. the circle of radius two) is marked so that one 
can clearly see the movement of the agents along this boundary. 
Figure El shows a top-down view of Fig. in which the pro- 
jected boundary of the surface (that is the circle of radius two) 
is marked and this enables one to see the accuracy of the motion 



12 



of the swarmers on this boundary. 



In Fig. [71 we pick one of the vehicles shown in Fig. [5] and 

show the relative error evolution for this particular trajectory. 
The chosen trajectory is the one that starts at xq = 0.4cos(^), 
Ho = 0.4sm(^), and it shows a relative error which is below 2%. 
The rest of the vehicles behave similarly. 

Developing the algorithm further, in Eq.( |Sj), we allow both the 
concentration C and the boundary Co to be functions of time as 
well as space. To simulate variation in time we consider the same 
surface as before given by ( ED and vary both x and y periodically 
in time with different frequencies according to: 

x = Aocos[^) y = Aosin{^). 

This results in a generic, time-dependent change of the whole 
surface. Here we show in Fig. |H1 four images of the motion at 
various moments in time. In Fig. IHl we show the relative error 
associated with the motion of one of the trajectories in Fig. El 
We chose again the same initial conditions as for Fig. [71 This 
error averages out to 3.6% when = 0.1 and stays below 5%. If 
we raise = 0.15, the corresponding average error will be 5.3%. 



5 Collective motion in 3D 



The above algorithm can be realized in 3D by a similar variational 
formulation. Initially the minimization problem was introduced 
in image processing [TS| as a deformable surface model meant to 
evolve to a certain location and shape. We extend that idea to 
the motion of a collective of agents. We consider the group of 
agents placed on a deformable surface which evolves in time to 
fit a certain boundary surface. The deformation assumes elastic 
and bending forces that are internal to the deformable surface. 
The forces have to balance the so-called external forces that are 
stretching this surface so that it will fit a boundary surface or 
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t=0.03 



t=0.23 





Fig. 5. The picture illustrates a swarm of 25 vehicles at different instances in time. 
The algorithm accurately follows a boundary variable in space. The motion is shown 
at t = 0.03 then at t = 0.23, 0.43 and t = 0.63. 



some other type of dynamical object that the moving swarms 
have to reach. We adopt, as in [T^, the model consisting of min- 
imizing an energy functional E over a set of parameterized ad- 
missible surfaces v defined by the mapping: 
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Fig. 6. The picture gives a top-down view of the swarm in Fig. 5. The algorithm 
accurately follows a boundary variable in space which, when projected forms a circle 
of radius 2. 



Q = [0, 1] X [0, 1] ^ 

(9) 

v{s, r) = {x{s, r), y{s, r), z{s, r)). 



The associated energy E is given by: 
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t 



Vk^_ ,^ . ^ ^ 

"0 100 200 300 400 500 600 

iterate 

Fig. 7. The picture shows the relative error corresponding to one of the trajectories 
in Fig. |SJ which is seen to stay below 2% 

E{v) = J^[wio\vs{s)f + woi\vr{s)f + 2wn\vsr{s)f 

+w^2o|^ss(s)|^ + wo2\vrr{s)\'^ + P{v{s , r))]dsdr (10) 

where P is associated to the task to be achieved by the algo- 
rithm. The boundary is now a 2D surface in space. Similar to 
the 2D case we let P{x,y,z) = {C{x,y,z) — Cq)^ where now C 
and Cq are functions of x,y and z. The other terms in Eq. [TUl 
determine how the motion in space of the surface v will take 
place. Continuing our analogy from continuum mechanics, these 
terms model internal forces that are acting on the shape of the 
surface v as it moves towards its goal. Consequently the shape 
of the surface v will depend on the elasticity coefficients wio, woi, 
the rigidity coefficients W2o^wo2 and the resistance to twist wn. 
We also constrain the surface v by assuming periodic boundary 
conditions. 

A local minimum v of E satisfies the associated Euler-Lagrange 
equation: 

-WiqVss{s, r) - WQiVrr{s, r) + 2wiiVssrr{s, v) + W20Vssss{s, r) 

+WQ2Vrrrr{s, t) = -VP{v{s, r)) (11) 

to which boundary conditions are added. Equation [TT] represents 
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Fig. 8. The picture shows four instances in time of the motion of a swarm of vehicles 
on a time-varying surface. The motion is shown at t = 0.3,0.7,0.11,0.15, the time 
step was At = 0.01 

the necessary condition for a minimum of {E'{v) = 0). A 
solution of Eq. [TT] can be seen as either reahzing the equihbrium 
between internal and external forces or as reaching a minimum of 
the energy E. Since there may be many local minima, appropriate 
initial data vo{s,r), leading to the desired minimum is assumed 
when solving the associated evolution equation, in which we add 
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iterate 

Fig. 9. The relative error in the motion of a single vehicle chosen at random out of 
the 25 trajectories in Fig. |HlThe error averages 3.6%. 

a temporal parameter t: 
dv 

— - WioVss{s) - WoiVrr{s) + 2wuVssrr{s) + W20Vssss{s) 

+W02Vrrrr{s) = -VP{v{s,r)). (12) 

A solution to the static minimization problem is obtained when 
the solution v converges as t tends to infinity, which means that 
the term ^ vanishes leading to a solution of the static problem. 

The evolution Eq. ( [T2j) will be used to model the motion of a 
swarm in 3D. We discretize this equation in space using finite- 
differences and each discrete spatial point will stand for the po- 
sition of a vehicle proceeding analogously to the 2D case in the 
previous section. Thus this swarm will be designed as lying on the 
surface v{s,r) and moving towards the assigned boundary with 
interactions between nearby vehicles determined by the elasticity 
and rigidity coefficients in Eq. (USD- We apply this discretization 
(not presented here, but standard and analogous to the discretiza- 
tion in section IV) to the motion of a specific swarm evolving on 
a concentration volume defined in 3D of equation: 
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C{x,y,z) 



tanh{F{x, z) 



3 



4 



(13) 



F{x,y,z) 



4 



{{x - x,)^ + {y - y,f + {z - z,)^) 

(t2 



) 



where {x,, y,) = (1, 0, -1), (0, -i i), (-|, i, 0), 1, -|) for i = 
1,2,3,4 and a = 1.0. 



The swarm will move in space through a medium of variable 
concentration in space. The swarm will move along the steepest 
descent direction which is given now by a 3-dimensional vector 

rdP dP dP\ 
^ dx^ i9y ' dz '' 

A last numerical example shows that the motion of the swarm 
in 3D can be directed to reach and track a specific curve on 
the boundary. We define a curve in space lying on the spherical 
boundary. This curve is defined by the equation: 

X = cos{(f)) cos{(t)) 

y = cos{(j)) sin{(j)) (14) 
z = sin{(t)) 

When using the discretization of Eq. [TT] with this modified goal 
dynamics we obtain the 14 trajectories seen in Fig. 10. The goal 
dynamics is introduced in the discretization of Eq. [TT] by evalu- 
ating Co at x^y and z given by [TH In Fig. 10 we show both the 
exact target curve (in red) a and the curve described by 14 swarm- 
ers (in blue) as they approximate movement along the curve of 
Eq. [HI In this figure we remark that the target curve is traveled 
by the swarmers in opposite directions. We also show a detail of 
Fig. [im in Fig. [TT] where arrows indicate how the trajectories 
evolve from inside the sphere towards the surface boundary and 
then move along the prescribed curve. 

As in the 2D case we now vary the concentration in time. In the 
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concentration equation, Eq. [T31 we let x, y and z vary in time 
according to the formula: 



x = Aqcos{^) y = AQCOs{^) z = Aqcos{^) (15) 

The corresponding figure is not shown here since it is very similar 
to Fig. 10 with the exception that the sphere undergoes slight 
deformations in time due to the motions given in [THl To illustrate 
the time-dependent case we give in Fig. 12 the relative error for 
one agent. Namely the one starting at: xq = 0.8cos(^)cos(8^), 
Uq = O.Scos{^)sin{8^), zq = 0.8sm(^). It can be seen in this 
figure that the relative error is below 5%. 



6 Conclusions 



We have presented an algorithm that allows a formation of agents 
to move towards a prescribed boundary and then describe that 
boundary. The motion takes place on stationary or non-stationary 
surfaces and can be realized as dynamics on a surface (2D) or as 
motion in space (3D). 

The algorithm represents a new development over snake algo- 
rithms used in image processing Those algorithms were 
primarily designed to detect level-sets in an image. Here we show 
that a similar algorithm will detect boundaries that vary in space 
and in time. Our new algorithm, along with the one in [i2j, di- 
rects the idea of snake algorithms towards solving a new problem, 
namely the problem of controlling motions of swarms. Another 
development was to allow the environment in which the swarm 
is moving to depend on time in a quite general manner. 

In future work we will extend the current algorithm to detecting 
curves with a more complex topology than the one considered 
here. Our main goal in future developments of these algorithms 
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Fig. 10. 14 trajectories starting inside the spherical domain, reaching the boundary 
and following a specified curve on the boundary. 




Fig. 11. Blowup of Fig. IIOI where the arrows show the direction of the agent's 
trajectories. 

will be to couple the equations of motion of the swarm with 
meaningful dynamical systems with relevance for biological ap- 
plication, chemical reactions and advective flows. In this way the 
abstract goal of reaching a certain surface or feature will be di- 
rected towards real-life applications. 
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Fig. 12. The relative error in the motion of a single swarm reaching and then moving 
along a space and time-dependent boundary in 3D. 
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